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Objective 

Currently no software simulations are available from industry or 
academia which can allow an accurate computation of the radar signature 
from full scale non-metallic vehicles primarily because existing 
formulations have excessive memory and computational demands in 
addition to being unsuitable for modern machine architectures. Such 
software are required for test and evaluation purposes of future aircraft 
systems/platforms, target simulation platforms, on-board communication 
systems and for radar signature evaluati on . 

Traditional methodologies for radar signature or performance 
analysis of on-board communication systems have not shown any promise 
in simulating large scale composite ve hi cles. Furthermore, approximate 
simulation techniques are of no utility in modeling complex, composite 
structures. The intent of this research is thus aimed at developing new 
methodologies which will allow accurate and efficient simulations of full 
scale vehicles. We propose to accomplish this with the development of 
algorithms which have very low memory requirements and make use of 
recent advances in high performance and parallel computing. The 
particular technique to be employed is a hybrid version of the finite 
element method which has been successfully implemented and tested at the 
University of Michigan for two-dimensional applications. Also, some three 
dimensional results relating to specific scattering and antenna 
configurations have verified our expectations of the method's potential. 
The basis of this approach is to reformulate an open-boundary problem 
into a closed-boundary one using a mesh termination method which brings 
the outer boundary very close to the body. Segmentation and connectivity 
schemes can then be employed to subdivide the system/volume into smaller 
(of any number) sections which can be independently treated by separate 
single or parallel processors. Except for the small Toeplitz connectivity 
matrices, all others are sparse and banded, making their solution 
convenient, and since there is no need for concurrent storage of any of 
these submatrices, there is no limitation on the size of the computational 
domain, a problem which has become the bottleneck in computational 
electromagnetics. It is thus our claim that this solution approach, coupled 
with the advantages of the finite element method, will make a dramatic 
impact in the field of numerical electromagnetics, and promises to provide 
the necessary future tools for numerous applications, including target 
simulation, radar signature analysis and imaging, antenna modeling and 
VLSI circuit simulations. 
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Progress 

As stated above, the main goal of this project is to develop methodologies 
for scattering by airborne composite vehicles. Although our primary focus 
continues to be the development of a general purpose code for analyzing 
the entire structure as a single unit, a number of other tasks are also 
pursued in parallel with this effort. These tasks, are important in testing 
the overall approach and in developing suitable models for materials 
coatings, junctions and, more generally, in assessing the effectiveness of the 
various parts comprising the final code. 

Below we briefly discuss our progress on the five different tasks which 
were pursued during this period. Our progress on each of these tasks is 
described in the detailed reports (listed at the end of this report) and the 
memoranda included with this document. The first task described below is, 
of course, the core of this project and deals with the development of the 
overall code. Undoubtedly, it is the outcome of the research which was 
funded by NASA-Ames and the Navy over the past three years. 


Task: Three-dimensional hybrid finite element formulation for 
scattering 

During this year we developed the first finite element code for scattering 
by structures of arbitrary shape and composition. The code employs a new 
absorbing boundary condition which allows termination of the finite 

element mesh only 0.3^. from the outer surface of the target. This leads to 
a remarkable reduction of the mesh size and is a unique feature of the code. 
Other unique features of this code include capabilities to model resistive 
sheets, impedance sheets and anisotropic materials. This last capability is 
the latest feature of the code and is still under development. 

The code has been extensively validated for a number of composite 
geometries and some examples are given in Figures 1 and 2. The 
validation of the code is still in progress for anisotropic and larger non- 
metallic geometries and cavities. 

The developed finite element code is based on a Galerkin's formulation and 
employs edge-based tetrahedral elements for discretizing the dielectric 
sections and the region between the target and the outer mesh termination 
boundary (ATB). This boundary is placed in conformity with the target's 
outer surface thus resulting in additional reduction of the unknown count. 
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The resulting system is symmetric and is solved via an iterative biconjugate 
gradient algorithm. It is particularly important to note that only the non- 
zero elements of the sparse matrix system. In this manner, the achieved 
storage requirement is 30 N. The edge based elements lead to matrices 
which have greater bandwidths but are sparser and the employed storage 
scheme takes advantage of this sparsity. In contrast to the node-based 
elements which have been traditionally used in aerodynamics, the edge- 
based elements are better suited for electromagnetic computations. These 
Witney form-1 elements maintain tangential field continuity and avoid 
specification of the field at metallic comers. This is crucial because certain 
field components become undefined at metallic comers. 

Over the last few months we have concentrated on the parallelization of the 
developed finite element code FEMATS. More specifically, the code was 
ported on our recently acquired Kendall Square multiprocessor facility 
equipped with 32 processors. At the moment, we have achieved a 23-fold 
speed-up demonstrating that the code is fully parallelizable using syntactic 
means. In the future, we shall couple syntactic parallelization with other 
means of parallelization based on the physical characteristics of the 
problem and associated formulation. 

In addition to the aforementioned unique analytical features of the code, we 
have also concentrated on user-oriented interface features. Specifically, the 
code is currently interfaced with the commercial mesh generation package 
IDEAS marketed by SDRC. This package has graphical interface for 
building and meshing geometries of arbitrary shape with look-up tables for 
material specifications. The Universal geometry file generated by IDEAS 
is read by a preprocessor which recasts the geometry data in a form 
suitable for the solver. Various postprocessors can be interfaced with the 
solver for displaying the near-zone fields and far-zone patterns in color 
coded form for diagnostic purposes. Postprocessors included in the IDEAS 
package, apE or grafic have been used for dat a display. To demonstrate 
this capability a 2-minute video animation was generated displaying the 
near and far zone fields of an aircraft wing as a function of incidence 
angle. This video was distributed to NAS A- Ames and Pacific Missile Test 
Center technical monitors of this project. 


Task: Hybrid finite element method for a body of revolution 

The body of revolution code was completed this year with the delivery of a 
detailed technical report 025921 -3 1-T describing the formulation and 
results based on the developed code. Several scattering patterns for coated 
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and uncoated bodies of revolutions are included in this report which 
validate the code and demonstrate its capabilities. A unique feature of this 
code is the combined finite element and boundary integral formulations 
employed in the implementation. The finite element formulations permit 
modeling of inhomogeneous materials, whereas the boundary integral acts 
as an exact mesh termination scheme. The mesh can therefore be 
terminated very close to the target thus minimizing the unknown count. 
The drawback of the boundary integral scheme is that it leads to a full 
submatrix which must be handled carefully for memory reductions. As 
done in our previous implementations of the finite element method, the 
boundary integral is cast in convolutional form and this leads to an 0{n) 
memory demand of the overall FE-BI system. Much effort was devoted to 
the careful computation and implementation of the boundary elements and 
details pertaining to this are described in the aforementioned report. More 
recently we concentrated on the application of the code to scattering by 
coated BOR structures and some examples are shown in figures 3 and 4. In 
particular, figure 3 displays the effectiveness of a recently developed 
scheme for suppressing the exterior/interior resonances associated with any 
solution employing a closed boundary integral formulation. Figure 4 
displays the scattering pattern by a coated cone frustrum. As seen the 
results of our FE-BI code are in good agreement with the moment method. 
A particular characteristic and advantage of the FE-BI method is its 
capability to easily treat surface details and inhomogeneous materials. 


Task: New impedance wedge diffraction coefficients 

Sometime ago we reported in the University of Michigan report 025921-4- 
T a new approximate impedance wedge diffraction coefficient for 
characterizing the diffraction by an impedance wedge. We resorted to an 
approximate coefficient only because the impedance wedge is not amenable 
to an exact solution when excited by a wave incident in a plane other than 
that normal to the wedge. This diffraction coefficient was later 
implemented and successfully tested in a traditional PTD code under a 
contract supported by the U.S. Army (see University of Michigan report 
028371-1-F). 

After several applications, the aforementioned diffraction coefficient was 
found to be at times inaccurate for computing the cross polarization fields. 
To remedy this deficiency, we examined the possibility of a more rigorous 
development for the proposed diffraction coefficients. From the start, this 
effort proved challenging but we have now arrived at a procedure which 
has yielded a more accurate diffraction coefficient. The new diffraction 
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coefficient is composed of three different ones, each valid at certain 
regions of interest as illustrated in figure 5. The validity regime of each of 
these is based on the specific approximations made in the derivation process 
and details pertaining to this as given in the technical report 025921-30-T. 
A chapter is also included in this report providing an extensive validation 
of the diffraction coefficients for different wedge angles and face 
impedances. These comparisons clearly show the derived diffraction 
tensors are very accurate. Examples of the comparisons are shown in 
figures 6 and 7. On these figures, die new solution is termed as “1st Order 
Approx. GTD”. The indicated differ e nces observed in the shown 
comparisons are attributed to multiple interactions, and do not therefore 
burden the developed diffraction coef ficien ts. It is worth pointing out that 
the derived diffraction coefficients reduce to the known tensors for the 
wedges having included angles of 0 and 180 degrees. 


Task: New finite element formulations for electromagnetic 

modeling 

In our last semi-annual report we reported of a finite element formulation 
using new solution-aware expansion functions. At that time, we had 
completed a two-dimensional implementation of this formulation. The 
three-dimensional implementation is now nearly complete and we shall 
intially consider its application to inlet and guided structures. For inlets 
having canonical cross-sections, the finite element method will be 
implemented with the guided modes serving as the basis functions in 
constructing the sparse system of equations. This should allow the 
simulation of long inlets since fine meshing will only be required at the 
inlet mouth, the engine face and at irregular bends of the inlet geometry 
not conforming to some canonical shape. This implementation will serve as 
a demonstration of the method's capability to model large inlet structures. 
The next step will be the implementation of the method to inlet structures 
whose shape may be considered as a perturbation of some canonical 
geometry. In the case of inlets having arbitrary cross section we will 
consider an implementation based on a combined finite element-high 
frequency implementation. That is, the fields in the guided region of the 
inlet structure will be formulated using a ray representation, whereas the 
regions near the inlet lip and the engine face will be, as before, fine meshed 
in the usual manner. 
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Task: New integral equations for scattering by coated 

structures 

Thin material coatings are nowadays routinely used in the construction of 
all airborne vehicles. Nevertheless, their electromagnetic characterization 
in the context of integral equations still relies on the use of the simple but 
less accurate impedance boundary conditions. On the other hand, although 
finite element formulations can readily deal with the presence of material 
coatings, their implementation is inefficient when the coating becomes thin. 
This is because very small elements must be used in the vicinity of the thin 
coating. 

To avoid the inaccuracy of the impedance boundary condition and the 
inefficiency of the finite element method for modeling thin coatings, a new 
integral equation was derived and discretized for implementation. The 
derivation of the new integral equation is described in the University of 
Michigan technical report 025921-25-T entitled “Alternative field 
representations and integral equations for modeling inhomogeneous 
dielectrics”. In comparison with traditional integral equations which 
require as many as eight unknowns in their discretization, the new integral 
equations requires only three unknowns per cell without any compromise 
in accuracy. The approximate integral equation based on the impedance 
boundary condition requires four unknowns per cell. In addition the new 
integral equation is also associated with a lower kernel singularity which 
greatly facilitates its discretization. A code based on this formulation has 
been written and validated. A short report pertaining to the formulation 
and results obtained from this code are described in the included 
memorandum. 
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( a/X 2 ) in dB 



Boundary 

Geometry of a cube (a=b=0.5A) consisting of a metallic section and 
an air-filled section, where the latter is bounded by a resistive surface 
having R = Z 0 . The outer cube drawn with dashed lines represents the 
surface on which the ABC was enforced. 



Figure 1. RCS patterns in the x-z plane for the shown geometry with H™ c = 0 
(solid curve is the FE-ABC pattern; circles are MOM data) and £‘ nc = 
0 (dashed curve is the FE-ABC pattern; black dots are MOM data). 
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(a/\*) in dB 



Observation Angle 6, deg. 


Figure 2: RCS pattern in the x-z plane for the geometry as in Figure 1, with the 
air-filled section replaced with dielectric having e r = 2 — j 2. The solid 
curve is the FE-ABC pattern and the black dots are MOM data for the 
E \ nc = 0 polarization. 
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Figure 3: The axial backscatter cross section is displayed as a function of normal- 

ized radius k 0 a o . The structure is a conducting sphere of radius a coated 
with a dielectric with e r = 2 - j4 has an outer radius of a 0 = 1.02. 
Employing a complex wavenumber k = k^ot , the resonance behavior is 
removed for the complex values in comparison to the a = 1 case 
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lap) 7 xi* 


The mesh of the sphere-capped cone frustrum 



(b) 


Figure 4: The bistatic scattering pattern is shown for a plane wave axially incid- 
ent on the coated sphere-capped cone frustum; H-plane pattern. 




Reflection SB. 




(b) 


Reflection 

S.B, 



Reflection 

S.B. 


Figure 5: Different regions for the applicability of the three solutions, (a) back- 
scatter case; (b) bistatic case. 
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1st Order Approx. GTD 
Moment Method 
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Figure 6: Far zone bistatic scattering pattern for a triangular cylinder with top 

side width w = 4A, a = 40°, & = 30°, 4>' = 70° , and its top face having 
t] = 0.5, and the other having tj = 1.0 + jl.0, (a) Ep 0 p' o , (b) £*<»'■ 
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Echo width/wavelength (dB) 


Figure 7: Far zone bistatic scattering pattern for a triangular cylinder with top 

side width w = 4A, a = 40°, (3 0 = 30°, <f>' = 70°, and its top face having 
tj = 0.5, and the other having r/ = 1.0 + j 1.0, (a) Eao> , (b) E^. 
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Operation count = 2.92 x 10 9 
MoM code 

No. of unknowns = 20,350 

Storage = 0.4 x 10 9 words 

Operation count = 8.4 x 10 12 



Number of operations vs surface area 
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Abstract 

A new integral equation for scattering by an inhomogeneous di- 
electric is discretized using rectangular brick elements. Edge-based 
Unear shape functions are used for the expansion of the field inside 
the dielectric and the modified Galerkin’s technique is employed for 
testing the integral equation. Results demonstrating the validity of 
the integral equation are presented. 


1 Introduction 

The modeling of inhomogeneous dielectrics by an integral equation approach 
is traditionally accomplished by the introduction of equivalent volume electric 
and magnetic currents as shown in figure (1). 
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Figure 1: Traditional Dielectric/Ferrite body formulation 

For a dielectric with non-trivial permittivity and permeability this type 
of modeling implies six scalar unknowns at each volume location. As a re- 
sult, the implementation of the resulting integral equation is computation- 
ally intensive and has excessive storage requirements. Recently [1] it was 
demonstrated that any inhomogeneous dielectric material, regardless of its 
permittivity and permeability profile can be modelled by a single electric 
or magnetic current density. Alternatively, either the electric or magnetic 
fields within the dielectric can be used as the unknown quantities. However, 
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since this reduced-unknown integral equation involves the derivatives of the 
unknown quantities, a higher basis function is required for discretizing the 
resulting integral equation. Hence the use of linear shape functions. 


2 Discretization 

2.1 The Integral equation 

The volume electric-field integral equation to be discretized is of the form : 

E‘(r) = E(r) - ^JVGo(ry) x I] • x E(r') 

+V' x [(e r — l)E(r')]| dv' (1) 

where I = xx + yy + zz, G 0 ( r, r') is the free-space Green’s function, E ! is the 
incident field and Vj is the dielectric volume. 

From the properties of an identity dyad, equation (1) becomes 

E’(r) = E(r )+[[[ VG 0 (r,r')xf— V'xE(r') 

J J JVj ^ fJ'r 

+V' x [(c r — l)E(r')]| dv' (2) 

Assuming the material parameters are constant in each unit cell of discretiza- 
tion and the domain of integration is limited to one such cell (say cell ‘q’) 
then equation (2) is : 

E‘(r) = E(r) + {<„ - -L J jj £ VG„(r, r') x [V' x E(r')] dv' (3) 

where /.i Tq and e Tq are the relative permeability and permittivity within the 
qth cell and V q is the volume of the qth cell. 

From elementary calculus, 

VG„(r,r') = -(A+^)G„(r,r')fi (4) 
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Hence equation (3) becomes: 

m - E(r)q fr ,--1}///,( A + I)6'„(r, r ') 

R x [V' x E(r')j dv' (5) 

2.2 Basis functions 

The dielectric layer is assumed to be thin (t <C A). It is discretized by 
use of rectangular brick elements whose height is equal to the thickness of 
the dielectric. For the expansion of the field inside the dielectric we use 
edge-based expansion functions borrowed from the finite element domain [2], 
Edge-based expansion functions avoid explicit specification of the fields at 
the corners and edges. The geometry of the brick element with a local co- 
ordinate system is shown in figure (2). The electric field can be expanded 
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as : 


N ( 4 
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Expanding equation (5) in terms of one such vector expansion function : 

E‘(r) = W;(r)-{e„-^}/// Vt (A + i)G„(r,r') 

R x [v' x W; i (r')] dt/ (7) 
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2.3 Testing procedure 

Galerkin’s procedure involves choosing the weighting (testing) functions to 
be the same as basis functions in the inner product of the method-of-moments 
solution as follows : 

^ = JIi w ^ dv = Ilij ^ dv («) 

where w and f are the weighting and basis functions respectively, and this 
integration is over the whole volume of the dielectric. Since the basis func- 
tions are subdomain basis functions, the integration limits of equation (8) 
actually can be reduced to V p , the volume of the cell p where the weighting 
function w is defined. Forming the above scalar product on equation (7): 


III w;,(r) ■ E'(r)<fo = Jfj v w;,(r) ■ WJ(r ')dv 




ftrq 


JJJv r w ^- JJLX^+id 0 ^ 

R x [V' x W”(r')] dv'dv 


(9) 


where 6(p—q ) indicates a contribution within the self-cell alone. Equation (9) 
gives rise to a system of equations which can be written as: 


{Z mn }{<f> n } = {V m } 


where typical elements are given by : 

VZ = ]JJ V W^(r) • E'(r)dv 


ypq 
* mn 


= 7/X W'(r)-W;(r')rf» 


ti(P~<l) - ) <r q 

L ^q y 

R x [v' x W; i (r')] dv' 


( 10 ) 
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For self-cells the kernel of the integrand in equation (3) has a singularity 
which needs to be carefully evaluated. For self-cells we re-write the corre- 
sponding impedance matrix element from equation (3) as : 


^mn 


Jjjjy * WJ(r ')) ■ J]f Vr W> (r) x VGoCr, r') dv dv’ 


(11) 


The symmetric property of the free-space Green’s function gives 


VG 0 (r,r') = — V'Gotry) 


( 12 ) 


Using equation (12) and interchanging the order of integration in equa- 
tion (11): 

*=. = [//x„ w ^ (r) ■ w ~ ( r ° H £(p - + 

//X (V' x W”(r')) • jfj v [Wf„(r) - W;,(r')] x VG 0 (r,r')dvdv' 

_/ __L 

jjj v (V x WJ(r')) ■ jw^(r') x V Jj^Go^ivj dv' (13) 


An analytical expression for the integration of the free-space Green’s function 
over an equivalent spherical volume is derived based on [3]. This is of the 
form : 



jY—le ,k ° So ( s } n ^° s °) . - Sq cos(fc 0 ' s o)) + sin(A.- 0 so)e ,koV 

KqSq { \ ko ) 

(i + ») - (£ + ;*>)} (14) 


where rj is the radius of a sphere equal in volume to the rectangular brick 
and so is the distance between the center of the sphere and the source point, 
r'. The name modified Galerkin’s method (MGM) comes from the fact that 
the testing cells are spheres while the source cells continue to be the original 


34 


rectangular bricks. Note that the MGM is not equivalent to using spherical 
cells because we continue to use rectangular bricks as the domain for the 

source integration. 

Also, 


v 7//,. g » < r ' r '> 


dv = 


[k 0 s 0 cos(A'oSo) - sin(A.- 0 so)] 


^o• s o eJ *' 0,, 


jt ] [A 0 .s 0 cos(A 0 .s 0 ) — sin(A 0 so)] 1 
^o• s o eJ^O,, j 

’ x{x - x c ) + y{y - Vc) + gz - z c ) 
■So 


(15) 


where (x c ,y c , z c ) are the co-ordinates of the cen ter of the sphere. Note that as 
■s 0 — > 0 all the terms in equations (14) and (15) converge to a finite quantity. 


2.4 Excitation 

In all the above cases, the excitation is a plane wave given by 

E*(r) = [(a • 0‘) S' + (o • 4>') 4>'] e -jk ' r (16) 

where a = S' cos a + <j>' sin a is the polarization vector given by 

k 1 = — A - o (sin 9' cos fix -f sin 0 l sin gy + cos O'z'j (17) 

and 9' and <f)' are the usual unit vectors in the spherical system and are 
associated with the angles 9' and (f>'. 


2.5 Scattering parameters 

The scattered field can be obtained from equation (2) and (4) as 


E' 


'(r)=//^VC„(r, r ')x 


^ -V'xE(r')+V'x[(e r - l)E(r')]}di/ (18) 

Mr J 


Making the far field approximation for the magnitude R sa r and the phase 
Rttr-r- r', 


E'M = ^ t {*, - Jg V X V' X E(r')l dv’ (19) 


-jk 0 r N 


_1 

d’rq 
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Then the RCS is given by : 


a = lim Am' 2 

T— >OG 


|E«(r )| 3 

|E'(r )| 2 


(20) 


3 Results and Validation 

The implementation of the presented numerical solution is lengthy as is usu- 
ally the case with most three-dimensional numerical solutions. The validation 
of the code was done with data obtained from codes used to generate the re- 
sults in [4]. Figure (3) displays the backscatter patterns for a 0.2A long, 0.2A 
wide and 0.025A thick dielectric plate. Figures (4) and (5) display backscatter 
patterns for 0.1 A long, 0.1 A wide and 0.025A thick and 0.2A long, 0.2A wide 
and 0.025 A thick dielectric/magnetic. plates respectively. Figure (6) shows 
bistatic scattering patterns for a 0.2A long, 0.2A wide and 0.025A thick PEC 
block simulated by a lossy dielectric. 


4 Conclusion 

In this work we discretized and implemented a new, reduced unknown inte- 
gral equation for scattering by a dielectric layer. Results demonstrating the 
validity of the integral equation were presented. The impedance matrix gen- 
erated was found to be asymmetric since the testing procedure destroys the 
reciprocity of the integral equation. Hence iterative schemes like the CGFFT 
for rectangular plates could not be employed; thereby limiting the size of the 
structures investigated. Modified forms of this integral equation which result 
in symmetric impedance matrices are being investigated. 
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oA [dB] 



Incidence angle, ©^degrees) 

Figure 4: Backscatter RCS for 1=0.1 A, w=0.lA and t=0.025A dielec- 

tric/magnetic plate (e r = 7.4 —j 1.11, fi r = 1.4 — jO. 672) 
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oA [dB] 



Figure 5: Backscatter RCS for 1=0. 2A, w=0.2A and t=0.025A dielec- 

tric/magnetic plate (e r = 7.4 — jl.ll,^ r = 1.4 — j0.672) 
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[ap] 7 X/° 


-15.0 



Incidence angle, G^degrees) 

Figure 6: Bistatic RCS for 1=0. 2A, w=0.2A and t=0.025A PEC simulated by 
lossy dielectric (e r = 1 — jr 1 00, fi r = 1) 
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